Simulations of pure and doped 
low-dimensional spin-1/2 gapped systems 



Nicolas Laflorencie and Didier Poilblanc 

Laboratoire de Physique Theorique, CNRS-UMR5152 Universite Paul Sabatier, 
F-31062 Toulouse, France 

didier . poilblancOirsamc . ups-tlse . f r , nicolas . laf lorencieOirsamc . ups-tlse . f r 

1 Introduction 

Many systems of Condensed Matter consist of fermions (electrons) moving on 
a lattice and experiencing strong repulsive interactions 1 . In such cases, the 
traditional perturbative methods to treat the electronic correlations often 
break down. In a pioneering work, Bonner and Fisher [2] revealed the ex- 
act diagonahsation (ED) method as a powerful tool to study the properties 
of one dimensional (ID) spin chains. Later on, it was extended to investi- 
gate two dimensional (2D) localised spin systems This work initiated a 
more extensive use of the method to investigate a wide variety of different 
systems such as strongly correlated lattice electrons (Hubbard-like models), 
mesoscopic systems 0], electron-phonon models, etc.... 

The success of this method first comes from the rapidly growing power of 
supercomputers which are being equipped with faster and faster processors 
and larger and larger memories, disk space and storage facilities. In addition, 
ED are clearly unbiased methods as we shall discuss later on in the course of 
this Chapter. The systematic errors can be, in most cases, easily estimated 
and, hence, this method is a very controlled one. Of course, it has its limi- 
tations (which we shall also discuss later) but, clearly, the efficiency of this 
technique will steadily increase in the future as the power of supercomputers 
booms up. 

The following Chapter will be dedicated mostly to the numerical technique 
based on the Lanczos algorithm. However, we shall focus on a specific area of 
strongly correlated models, namely low-dimensional spin-^ gapped systems, 
to illustrate various technical aspects of the method and to discuss the physics 
of these topics. References to related specialised work dealing in more details 
with the physical aspects will also be given. 

Low dimensional spin-^ systems with antiferromagnetic (AF) interactions 
display very innovative features, driven by strong quantum fluctuations. In 
particular, geometrical effects or competing magnetic interactions can give 
rise to the formation of a spin gap between the singlet ground state and the 
first excited triplet state. In this chapter, we focus on the numerical inves- 
tigation of such systems by Exact Diagonahsation (ED) methods and some 
extensions of it including a simultaneous mean-field (MF) treatment of some 
perturbative couplings. 
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This Chapter is organised as follows: in Section 2 a description of the Lanczos 
algorithm is given with special emphasis on the practical use of space group 
symmetries. A very short review on the well-known planar frustrated Heisen- 
berg model and some linear chain Heisenberg models is given in Section 3. In 
particular, we outline the role of the magnetic frustration in the formation 
of a disordered phase. We also introduce a MF treatment of interchain cou- 
plings. Section 4 is devoted to more recent studies of impurity doping and to 
the derivation of effective models describing interaction between dopants. 

2 Lanczos Algorithm 
2.1 Algorithm 

The exact diagonalisation method is based on the Lanczos algorithm 5 which 
we shall describe here. This algorithm is particularly suited to handle sparse 
matrices and there is, in fact, a wide variety of lattice models belonging to 
this class as we shall see later on IQ. 

Let us consider some lattice model corresponding to some Hamiltonian H 
with its symmetry group Q = {g}, namely [H^g] = 0. Let us also assume, for 
the moment, that irreducible representations of the symmetry group can be 
constructed. They consist of complete subsets of states Ai = {la)} which are 
globally invariant under the application of the Hamiltonian H (we postpone 
to the next part of this section the explicit construction of these states). 
Clearly, H can be diagonalized in each of the subsets Ai independently. It 
can be written as a tridiagonal matrix in a new orthonormal basis set {|??m)} 
defined as j2j 

H\<Pi) = ei\<Pi) + b2\<p2), 

■■ (1) 

H\<Pn) = e„|<?„) -|-6„+i|'^„+i) +bn\<Pn-l) , 

where the various coefficients and new basis states can be calculated recur- 
sively. The proof is as follows: let us suppose that the procedure has been 
applied until the step n, i.e. an orthonormal set of states |^i), \'Pn) has 
been constructed. Assuming the knowledge of ei,...,e„_i, &2,---,6n, |^n-i) 
and l^n)-, one can then determine 

e„ = i^nim^n). (2) 

Hence, the new state defined by 

\(j)n+l) = H\<Pn) ~ e„|<P„) - hn\^n-l), (3) 

is clearly orthogonal to |<?„). Moreover, (^„_i|0„+i) = (^„_i|iJ|^„) - 6„ 
which is also vanishing as can be seen by substituting the expression for 
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More generally, \(j>n+i) is, in fact, orthogonal to all the previous 
states \^p), p < n as can be shown recursively. Indeed, let us assume that, 
for p < n, 

Vi, p<i<n {<P,\(f>n+i) ^ 0, (4) 

then (<Z>p_i|0„+i) = {$p^i\H\<Pn) where (<Z>p_i|<?„) = {$p^i\^n-i) = has 
been used. Substituting the expression given by Eq. ^ for H\(pp-i) leads to 
the expected result 

(<?p.i|0„+i) =0. (5) 
The (positive) number 6„+i is simply defined as a normalisation factor, 

bl+l = (0n-|-l|0n+l), (6) 

i.e. |<?„+l) = ^J(j)n+l). 

In principal, a zero vector will be generated after iterating the Hamilto- 
nian a number of times corresponding to the size J\fi of the Hilbert space. 
However, the number of iterations necessary to obtain the lowest eigenval- 
ues and eigenvectors is much smaller. Typically, for A/j 10^, the ground 
state can be obtained with an accuracy better than 10^^, by truncating the 
procedure after only Na ^ 100 iterations and by diagonalizing the resulting 
tri-diagonal matrix by using a standard library subroutine. However, for a 
given size Afi of the Hilbert space, the necessary number Nu of iterations 
might vary by a factor of 2 or 3 depending on the model Hamiltonian. In 
practice, the convergence is faster for models for which high energy config- 
urations have been integrated out (e.g. t-J models in contrast to Hubbard 
models). Note however that, in some cases (models with strong finite range 
interaction) , the energy vs Na curve can exhibit steps before convergence to 
the true ground state is achieved. Once space group symmetries have been 
implemented, the best choice for the initial state is a purely random 
vector. The ground state is also easily obtained as a function of the states 
|<?ri). However, to express it in terms of the initial basis {\ct)} (as it is often 
useful) it becomes necessary to store temporarily the intermediate vectors 
|<?„). This step is usually the most demanding in terms of disk space and/or 
mass storage. However, note that, at runtime, only three Lanczos vectors of 
size J\fi need to be assigned in memory. 

Full diaqonalisation: In some special cases (spectrum statistics thermo- 
dynamics j9,, etc ...) the complete spectrum (or a least the lower part of 
it) is needed. This can also be performed by the Lanczos algorithm. In this 
case, usage of a more sophisticated standard library package (e.g. EA15 of 
Harwell) is preferable. Indeed, more involved tests become then necessary to 
eliminate the unphysical "ghost" levels appearing (always above the ground 
state energy) after the diagonalisation of the tridiagonal matrix. However, 
the input for the library subroutine consists only on the set of states 
which have to be calculated separately (see below). 
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2.2 Space group symmetries 

Usually, some efforts have to be carried out in order to take full advantage 
of the symmetries of the Hamiltonian. Let us consider a model defined on 
a JD-dimensional lattice describing interacting fermions as e.g. the simple 
Heisenberg (spin) model, 

Hj = ^ Jxy Sx • Sy, (7) 

where interaction (in this case the exchange coupling) is not necessarily re- 
stricted to nearest neighbor (NN) sites but, nevertheless, is supposed to ex- 
hibit translation and point group symmetry. In other words, denoting the 
point group^ by Qp = {.gp}, we assume, e.g. in the case of Eq. ((TJ, 

Jxy ^(x-y) 

and (8) 
VffP e Gp, J{gp{v))^J{v). 

Clearly, such properties are easily generalised to generic spin or fermionic 
Hamiltonians. In addition, we shall also assume, as in Eq. (|7|), spin rotational 
invariance (the total spin S* is a good quantum number) or, at least, in- 
variance of the Hamiltonian under a spin rotation around some quantisation 
axis. Translation symmetry can be preserved on finite systems provided the 
geometry is that of a £>-dimensional torus with periodic or twisted boundary 
conditions (BC). On the torus geometry, the full space group reads, 

g = gp<E)T, (9) 

where we denote by T = {tp}, p = 1, TV, the translation group. It is clear 
that Hj is invariant under any g (z Q. 




Fig. 1. Schematic representation of of different lattices, (a) the y/32 x square 
lattice with NN and second NN couplings Ji and J2. (b) the dimerized ring with 
dimerization S. (c) The frustrated ring with NN and second NN couplings ,7i and 
J2. 
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2.3 Construction of the Hilbert space 

The motivation to take into account the Hamiltonian symmetries is two-fold. 
First, the Hamiltonian can be block diagonalized, each block corresponding 
to an irreducible representation of the symmetry group. Practically, the sizes 
of the various blocks Ni are much smaller than the size of the full Hilbert 
space (see e.g. Table^, hence, minimising the numerical effort to diagonalize 
the Hamiltonian matrix. Secondly, each irreducible representation of the sym- 
metry group is characterised by quantum numbers (such as the momentum 
k) which can be connected directly to physical properties. 

Table 1. Symmetry groups and sizes of the reduced Hilbert spaces for various spin- 
i AF Heisenberg models for one of the most symmetric irreducible representation 
(typical reduced size is written as the total size in the = sector divided by 
the number of symmetries). The ID models are descibed in section [3.21 Tjv and 
I2 stand for the translation group T (see text) and the spin inversion symmetry 
S'x — + — , respectively. 



Model 


Lattice Size 


Symmetry group 


Typical reduced Hilbert space size 


2D Isotropic 
ID Ji - Ja 

ID Ji -J2-5 


6x6 
32 X 1 
32 X 1 


Tae ® Civ ® h 
T32 ®C2®h 
Tie ® h 


9 075 135 300/576 
601 080 390/128 
601 080 390/32 



We shall focus here on the Heisenberg model but t-J and Hubbard 
models can also be studied. In order to minimise memory occupation, con- 
figurations can be stored in binary form. One, first, assumes an arbitrary 
labelling of the lattice sites from 1 to L and denotes a configuration in real 
space by 

|c) = |S1, ...,Si, ...,sl) . (10) 

For the Heisenberg model, the information Si on each lattice site can be 
stored on a single bit, e.g. "f and | correspond to 1 and respectively. A 
spin configuration with up to 64 spins can then be represented by a single 
64-bit word. As an example, a = 4 site configuration Heisenberg model 
such as I t, i, T, i) is coded as O64...O5I4O3I2O1 (where the subscripts indicate 
the place of the bits) i.e. the integer "10". An integer N{\c)) can then be 
uniquely associated to each configuration |c) and formally written as, 

N 

Ar(|c)) = ^2-V. , (11) 

1 

for spin-1/2 models (ct^ = 0, or 1). 

At this point, it becomes useful to consider space group symmetry. Each 
irreducible representation can be characterised by a momentum 
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K = ^n^K^, (12) 
ft 

where are the reciprocal lattice vectors (e.g., in 2D, = ^T^^, A 62) 
and are integers. For each value of K, one then considers G^, the so-called 
little group of K {Q^ <Z Gp), containing group elements gp such that 

5P(K)=K. (13) 

The relevant subgroup of Q to be considered is then 

Gk = GS:®T. (14) 

For a given symmetry sector / = (K,tk) (tk is one of the irreducible rep- 
resentations of Gk) * "symmetric" state |a) can then be constructed from a 
single configuration |c) as a linear combination which reads, up to a normal- 
isation factor, 

\a) = \a){\c)}= e(rK,5p)exp(iK-Tt) (5pi)(|c)), (15) 

gpegP.teT 

where e{TK,gp) are the characters (tabulated) of the representation tk and 
Tt are the translation vectors associated to the translations t. Since the 
procedure to construct the symmetric state is well defined, it is clear that 
one needs to keep only a single one of the related configurations gpt{\c)), this 
state being called "representative" of the symmetric state \a) and denoted 
by |r) = R{\c)). This naturally implies, 

V^e^K, R{g{\c))) = \r). (16) 

In other words, one retains only the configurations |c) which can not be 
related to each other by any space symmetry g G Gk- The set of all the rep- 
resentatives (labelled from 1 to A/;) defines then unambiguously the Hilbert 
space Ai = {|a)}. Typically, the size of this symmetric subspace is reduced 
by a factor of card(5K) compared to the size of the original basis |c). Note 
that, in some cases, there might exist some configurations |c) (to be elimi- 
nated) which do not give rise to any representative. This occurs when there 
is a subset Q' C Gk such that 

e G', g{\c)) = \c) 

and (17) 

^ e(rK,5p(5))exp(iK-T(5)) =0. 

g&S' 

The choice of the representative among the set of equivalent states (i.e. 
states related by some symmetry operations of Gk.) is, in principle, arbitrary. 
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However, as we shall see in the following, it is convenient to define it as the 
state |c) of a given class giving rise to the smallest integer N{\c)) i.e. 

iV(|r))= min{iV(g(|c)))} (18) 

For simplicity, we shall here extend our coding procedure to more general 
Hubbard-like models where one can construct the configurations by a tenso- 
rial product of the up and down spin parts 

|c> = |c(T)> ® |c(i)) . (19) 

N{\c)) contains now 2N bits and is of the form 

7V(|c)) = 7V'(|c(T)))x2^-|-iV'(|c(i))), (20) 

where N'{\c{(t))) corresponds to the binary coding {N' {\c{a))) < 2^)) of the 
cr-spin part of the configuration, by using Eq. (|11() where 1 (resp. 0) refers 
now to the presence (resp. absence) of a spin (7(=t oil) at a given site i. 
Although this coding is more costly in term of memory space (2N bits per 
configuration instead of N) , it is more general and applies both to Hubbard- 
like models (where "f and | spins can leave on the same site) and to Heisenberg 
models. The minimisation of N{g{\c))) over g E t/K can be done in two steps. 
First, one generates all possible configurations for the up spins (this usually 
involves a limited number of states) and only representatives jj'd)) are kept. 
At this stage, one needs to keep track of the subsets of symmetries fKi^^d))] 
of Qk leaving these representatives invariant. In a second step, one constructs 
the full set of configurations as a tensorial product of the form \r{1)) (g) |c(4)). 
The remaining symmetries of [!''(?))] ^'^e then applied to the spin J, part 
and one only retains the configurations |c(|)) such that, 

Vg' e £K[r(T))] N'i\ca))) < N'[g\\cim . (21) 

The Hilbert space is then formally defined by all the symmetric states la)!!?")} 
(see Ea.()15|lV However, one only needs to store the binary codes N{\r)) as 
well as the normalisations of the related symmetric states. The normalisation 
of H15|l requires some caution. In some cases, there might exist more than a 
single (i.e the identity operator I) group element g' of £k[|?'(T))] which keeps 
|r) invariant. Then, 

.FK[|r)] = {g' e £K(|r(T)));iV'[5'(|c(i)))] = A^'dcQ)))} (22) 

defines the subgroup of such symmetry operations. The sum over the group 
elements in Eq. H15|) should, in fact, be restricted to the elements of the coset 
Gk/J'-kIIt)] and the appropriate normalisation factor is then given by 



(23) 
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The two integers N{\r)) and n{\r)) corresponding to the binary code of a 
representative and to the normahsation of the related symmetric state, re- 
spectively, can be combined and stored in a single 64-bit computer word. 

In the case of the generic Heisenberg model, the previous two-step proce- 
dure can, in most cases, also be done using the N bits coding from Ea. (|ll|l 
but it is more involved. We only indicate here the spirit of the method. More 
technical details for spin-1/2 models can be found in Refs. JT] and The 
decomposition between up and down spins is replaced here by an (appropri- 
ate) partition of the lattice sites into two subsets A and B. The computer 
word (integer) coding each configuration contains then two parts, each 
part corresponding to one subset of the lattice sites, \c{A)) and \c{B)). The 
partition is chosen in a way such that the group can be decomposed as 

aK = a^+5(g|), (24) 

where 5 is a group symmetry which fulfils = I, the subgroup Gk leaves 
the two subsets of lattice sites globally invariant and S{G^) is the coset of S 
relative to G^- The decomposition (|24|) is not always unique. For 2D clusters 
such as the ^/32x ^/32 lattice of Fig.Q] a convenient choice for 5 is a reflection 
symmetry. Note that, in the case of ID rings, (|24|l is only possible if the 
number of sites is of the form = 4p -I- 2 or for very special values of the 
momentum K (like K — or K — n). Then, the previous procedure can be 
extended to this case by writing the configurations as 

|c) = \c{A)) ® \c(B)) (25) 

and by (i) applying the subgroup G^ on |c(^)) to generate its corresponding 
representative |f (A)) and then (ii) applying all the symmetries of £K[k(^))] 
on the part \c{B)). The action of the remaining symmetry S is considered at 
last; if the lattice sites are labelled in such a way that 

Xi+jv/2 5(xi) , (26) 

for i < N/2, the application of S can be implemented as a simple permutation 
of the two sub- words of N{\c)). 

2.4 Construction of the Hamiltonian matrix 

We turn now to the implementation of the basic operation |<25„) H\(l>n) 
appearing in Eq. ^ which is always specific to the model Hamiltonian H 
and constitutes the central part of the Lanczos code. Since the states \<Pn) are 
expressed in terms of the symmetric states \a) of (|15|l the Hamiltonian matrix 
has to be computed in this basis. For this purpose, it is only necessary to 
apply H on the set of representatives \r^) (labelled from 1 to A/;). In general, 
each configuration |r) leads to a small number (imax (at most equal to the 
number of terms in H) of generated states. 
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H\r,) JY.{-lt-^'^\c,,p) , (27) 

15=1 

where different signs (— l)^^^''' might arise (in the case of fermion models) 
from fermionic commutation relations. The matrix is then very sparse. Note 
that the full Hamiltonian can always be split in a small number of separate 
terms so that the amplitude of the matrix elements in H27|l is just a constant 
and hence does not need to be stored. 

At this point, it becomes necessary to determine the representatives (in 
binary form) of the various generated states on the right hand side of H27|l 
by applying to them all the symmetries of CJk- To achieve this, the choice of 
Eq. (|20|l for the binary coding of the configurations is very convenient. It is, 
indeed, a simple way to take advantage of the natural decomposition of the 
generated states, 

|c^,/3) = |c^,/5(T)) ® \c^Al)) ■ (28) 

Although, we restrict ourselves here, for sake of simplicity, to the general 
coding of the Hubbard-like models, the following procedure can be straight- 
forwardly applied using the more restrictive Heisenberg form H25|l for the 
configurations |c^,/3). The calculation of the representative 

\r^,p)f^R{\c^An®\c-^Al))} (29) 

can be done in two steps. First, one applies all the symmetries of to 
1^7, /3(T))- Since this procedure has to be repeated a large number of times, 
the function 

R : |c(T)) ^ |r(T)) . (30) 

can be, in fact, tabulated, prior to the actual calculation of the matrix el- 
ements. This is made possible since the number of possible states \c{])) re- 
mains, in general quite modest. This procedure enormously speeds up the 
calculation of the representatives and justifies the choice of Eq. H2()|l . Note 
that one also needs, in this preliminary calculation, to store, for each config- 
uration |c(t)), the corresponding ensembles, 

7^K[(|c(T))] = {.9 e ^k;.9(|c(T))) = k(T))} . (31) 

This also requires limited storage since, in most cases, 7?.K[|c(t))] contains 
just a single element. In a second step, it is sufficient to apply only the re- 
maining symmetries o^^Z\^^c^A'\))] to the |c-y^/3(|)) part. A standard hashing 
table is then used to find the positions of the representatives in the list. 
The connectivity matrix connecting the labels of the initial set of represen- 
tatives to the labels of the new set of generated states is then stored on 
disk. 

Note that the phases related to commutations of fermion operators and/or 
to the characters of the symmetry operations involved in the transformation 
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of the generated states to their representatives have also to be kept. These 
phases have the general form, 

= (-l)'-"'^e(TK, gp(5;^)) exp (zK • T(5;^)) , (32) 

where 5* ^ is a group element of 7?.K[|c-Y,/3(t))] such that 

k7,/3)/ = gjAl^-r.fi)) > (33) 

which depends on 7 and (3. It is easy to show that, if there exists more than 
a single group element which fulfils then, all of them lead to the same 
phase factor (|^ . 

Since the number of possible different phases given by (|32l) is quite small, 
it is possible to store the A^^^ (in some convenient integer form) together with 
■^{k7,/3) /} ^ small number of computer bits. For an Hilbert space of 10* 
(hundred millions) representatives with typically an average of ~ 50 images 
per state, the occupation of the disk is of the order of 5000 Mw i.e. 40 Gb. 
This can even be reduced by a factor of two by using each computer word to 
store the informations corresponding to two images instead of a single one. 
Once it has been generated, the Hamiltonian matrix is cut into several pieces 
(typically of the order of 10 to 100 Mw) and the various parts are successively 
read from the disk in order to calculate H\<Pn)- The best performances are 
obtained when the calculation using the n*'' part of the matrix and the access 
to the disk to read the (n+1)*'* part are simultaneous. Note that the Lanczos 
algorithm as it has been described above is well adapted to be implemented 
on a vector supercomputer (e.g. on a NEC-5 x 5). 

3 Examples of translationally invariant spin gapped 
systems 

Here, we first restrict ourselves to systems where the symmetry analysis de- 
cribed above can be used. Note however that explicit symmetry breaking may 
be present but, in general, the remaining symmetry group contains a large 
number of symmetries which can be exploited. Note also that since sponta- 
neous symmetry breaking can occur only in the thermodynamic limit, it does 
not prevent for finite systems the previous symmetry analysis. 

3.1 Application to the 2D Ji — J2 model 

The study of quantum phase transitions in 2D is of great interest. One of the 
standard example is the so-called frustrated Heisenberg model (see Fig^a)). 
It is defined by O when only NN and second NN exchange couplings are 
retained, 
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J^y = Ji if |x - y| 1, 

Jxy - J2 if |x - y| = \/2, (34) 
J^y = otherwise . 

Various analytical approaches have been applied to this problem such as spin- 
wave calculations large-N SU(N) theories jJHli series expansions J^Ij and 
Schwinger boson mean field approaches J7]. However, most of these meth- 
ods are somewhat biased since they assume the existence of some particular 
ground states. Pioneering unbiased ED studies ^lElED] have strongly sug- 
gested the existence of a disordered magnetic phase for intermediate couplings 
J2/ J^- Here, we briefly discuss some more recent subsequent work |21II11[[T^ 
which attempted to obtain more accurate results by a finite size scaling analy- 
sis. Similar studies have also been performed for other S = 1/2 2D spin mod- 
els like the triangular lattice [221 7 the Kagome lattice [53], the 1/5-depleted 
square lattice or the (2D) pyrochlore lattice |25|. 




Fig. 2. Finite size results for M^{tv, tt) for different values of J2. The dashed lines 
are least squares fits to the data, using all available clusters. The full lines are fits 
using only A'' = 20, 32, 36 (Reprinted from Ref. )12,'l. 



A magnetic ordered phase with a spin modulation Q can be characterised 
by an order parameter Mn{Q) defined by 

1 ^ 

MUQ) = j;^7^^A%\iY.expit^, ■ Q)S.,f\%) , (35) 



1=1 
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where jtf'o) is the ground state. It is important to notice that, in any finite 
system, the order parameter itself has a zero expectation value in the ground 
state due to spin SU(2) symmetry. In other words, the macroscopic mag- 
netisation can slowly fluctuate so that in average it vanishes. It is therefore 
essential to consider, as in (|35|) . the square of the order parameter which can 
be interpreted as a generalised susceptibility. 

For weak frustration J2/ Ji, Neel order with Q = (tt, tt) is expected, while 
for large ratio J2/ Ji a collinear phase with Q = (tt, 0) or Q = (0, tt) consisting 
of successive alternating rows of parallel spins is a serious candidate. Indeed, 
in such a collinear phase, each sub-lattice has Neel order so it is clear that 
it is stabilised by J2. Note that the normalisation factor of the staggered 
magnetisations (13511 is chosen so that the order parameter is independent of 
the size in a perfect classical Neel or collinear state. In such ordered phases 
where the continuous spin symmetry is spontaneously broken, field theory 
arguments [201 suggest a scaling of the form, 

M^(Q) mo(Q) + ^ . (36) 




Fig. 3. Comparison of the finite size fits for the anti-ferromagnetic and collinear 
order parameters (left and right curves, respectively) with linear spin wave theory 
(Reprinted from Ref. [H]). 



Square clusters with N ^ Ap sites (so that (tt, 0) belongs to the reciprocal 
lattice) are considered (see Sec. II. B), i.e. N=16, 20, 32 and 36. As seen in 
Fig. 121 corresponding to Q = (7r,7r) the scaling law H36|l is very well satisfied. 
Note that the 4x4 cluster shows systematic deviations. Similar results are 
also obtained for the collinear order parameter at larger J2/J1 ratios. 
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The extrapolated results are shown in Fig. Q . The most interesting fea- 
ture is the existence of a narrow range of Jij J\ around 0.5 where none of 
the ordered states is stable. Various candidates for this disordered phase 
have been proposed such as the dimer phase |15j and investigated numeri- 
cally UHlEniEIl- More details on this topic can be found e.g. in the review 
article [TT] . 

3.2 Application to spin-Peierls chains 

a) Purely ID models 

We now move to systems with anisotropic couplings in space (quasi ID mate- 
rials). Let us first consider purely ID models in order to describe the physical 
origin of the spin-Peierls (SP) transition. One of the simplest ID model which 
diplays such a behavior is the frustrated spin-i ring (see Fig^c)), also called 
Ji — J2 or zig-zag chain, described by the Hamiltonian 

L 

Hhust — ^^(^iS^i-S'i+l + J2Si.Si+2)- (37) 

1=1 

Its symmetry properties are given in Tabled for L = 32 sites. 



0.246 
0.245 
^ 0.244 
^" 0.243 
0.242 
0.241 





^0 























0.000 0.005 0.010 0.015 0.020 
1/L^ 



Fig. 4. Critical value of the frustration vs the inverse square of the system size 
obtained using the method developped in |27| . (Reprinted from Ref.|28p. 

The low energy properties of such a model are very interesting because 
it is gapless as long as a — J2/J1 remains smaller than a critical value 
Qc — 0.2412 For a > ckc, a gap As{a) cx e~^"~"=) ^ develops and 
a spontaneous dimerization appears, characteristic of the SP transition. At 
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a — 0.5, the so-called Majumdar-Ghosh (MG) point j30 , the 2-fold degener- 
ate ground state is known exactly and consists in the product of spin singlets 
located either on odd or on even bonds. Beyond the MG point, the short- 
range correlations become incommensurate. The triplet (5 = 1) spectrum of 
the SP phase is a two-particle (so-called kink or soliton) continuum as evi- 
denced by the scaling of the soliton- antisoliton binding energy to zero |31| . 

Adding an explicit dimerization , i.e. a rigid modulation 6 of the NN 
coupling (see FigQ]), drives immediately the ground state into a SP phase for 
any 6^0 even if J2 = 0. In term of symmetries, C2 and the translations 
of an odd number of lattice spacings are lost (see Fig^b) and Tabled])). If 
J2 ^ 0, the dimerized and frustrated model 

L 

Hdnu = I] ^1 [(1 + i-iyS)S,.S.,+, + aS,.S,+2] (38) 

displays an enhancement of the dimerized gapped phase, indeed As{a,S) — 
Z\s(a,0) cx 52/3 



6 



Incommensurate 



SP 

Commensurate 



U gapless 



Fig. 5. Phase diagram of the frustrated dimerized Heisenbeg AF spin-i chain in 
the a — S plane. The dotted line is the Shastry-Sutherland line {2a + S = 1) |33j . 
On its left side, the phase is SP gapped and commensurate whereas on the right 
side, the correlations are incommensurate. 



These properties are summarized in the phase diagram shown in Fig|5| 
Note that the static modulation S leads to soliton- antisoliton boundstaes as 
shown in FigEl 

b ) Chain mean field theory for coupled spin chains 

Physically, the previous models are often inadequate to describe the proper- 
ties of several compounds like CuGeOaP^ or LiV205|35| which are excellent 
realizations of weakly interacting frustrated spin-i chains. Let us first con- 
sider a set of frustrated spin chains which are coupled by a weak AF exchange 
Jx. This 2D model is governed by the following Hamiltonian 
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1.2 

E_ 

J 

0.8 



0.4 



0.0 

" k 

Fig. 6. Lowest lying triplet(»), singlet (o) and quintuplet ((g)) excitations vs. the 
wave vector k for the dimerized frustrated chain Eg. 1381 with J2 = 0.5, S = 
0.05, L = 28. Results for L = 32, fc = are shown to the left. ED results reprinted 
from |36|. 

L M 
i—1 a—1 

+ J±Si^a-Si^a+l]- (39) 

where i is the lattice index along the chains of lenght L and a labels the M 
chains {L and M are chosen to be even and periodic boundary conditions are 
assumed in both directions) . Obviously it should be possible to study exactly 
this spin model on the square lattice but, as we have seen above, it is hard 
to perform ED with system larger than 36 spins. Here, we take advantage of 
the fact that J_l << 1 to perform a MF treatment of the transverse coupling. 
Following Schulz Tf, the chain mean-field (CMF) version of (I39|l is given by 




L M 
i—1 a—1 

+KaSU- Ja.{sIJ{si,^,)1 (40) 

with 

Ka = MiSla+l) + (SL-l)), (41) 

the local magnetic field to be computed self-consistently. In the absence of 
dopant (see Section 4), we expect an homogeneous AF phase characterized by 
a self-consistent staggered magnetization (S'f q) = (— l)*+°m. Therefore the 
coupled chains problem is reduced to a single chain in a staggered magnetic 
field = ±2{-iyjj_m. 
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L 

-ffsingio(a, Jl^ = ^[S',.S',+i + + 2mJj_{-iyS^] + constant, (42) 

1=1 

and the symmetry group of such a model is Tj^ /2- In the absence of frustration 
{a — 0), it was shown that m ^ VJi 37 . By solving the self-consistency 
condition using ED of finite chains the transition line Jj_ = J±{<x) (see 
Fig. [71) separating the dimerised SP phase (m — 0) and the AF ordered 
phase (for which m 7^ 0) has been obtained in agreement with field theoretic 
approaches '38'. Finite size effects are small in the gapped regime and espe- 
cially at the MG point. Note also that numerical data suggest that the AF 
order sets up at arbitrary small coupling when a < ac with a clear finite size 
scaling J1{L) cx l/L at small a. 




0.1 0.2 ^ 0.3 0.4 0.5 0.6 0.7 



a 

Fig. 7. Phase diagram of the coupled frustrated chains as a function of the frus- 
tration a and the inter-chain magnetic coupling J±. The points, calculated for 2 
different chain sizes (12 and 16 sites), separate a dimerized phase (SP) from a Neel 
ordered phase (AF). The closed diamond shows the order-disorder critical point at 
Qc — 0.2412. The dashed line represents the expected behavior in the thermody- 
namic limit. In the inset, the staggered magnetization m{J±) has been calculated 
for a L = 12 sites chain along the MG line (dot-dashed line). Along the a = Q line, 
different symbols show the critical J± for L = 8, 12, 16, 18 from top to bottom and 
we have checked its scaling to according to a l/L law. 



c) Convergence issues of the numerical CMF 

The numerical procedure consists of successive Lanczos-diagonalizations of 
a frustrated (a) spin-1/2 chain. At each step the AF order parameter m is 
calculated and reinjected at the following step as the "new field". Starting 
the numerical procedure with an arbitrary value of m(0) 7^ 0, the chain of 
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size L is first diagonalized, m(l) is extracted and then used for the next 
iteration. Eventually the procedure converges to the fixed point m* . A very 
interesting feature is that the convergence to m* as function of the number 
of MF iterations p (see FiglHl is exponential. 

m(p) — TO* cx exp(— p/^T-) forp>>^T-, (43) 

with S,t{J±) a typical convergence time scale. 

In order to study convergence at large p we have considered here small 
systems (12 sites) at the MG point where the finite size effects are very 
small. We have checked this convergence issue by studing the speed V{p) = 




Fig. 8. Universal behavior of the convergence speed of the MF iterative procedure 
plotted versus the renormalized iteration index t = j-. Results are shown for a 
system of spins interacting via Eg. 11421 with L = 12 and a = 0.5. From top to 
bottom J_i = 0.075, 0.13, 0.1, 0.11, 0.108, 0.106. An initial value m(0) = 1/2 is used 
for all simulations. Convergences to m* = if the phase is SP (solid lines) or m* 7^ 
(long-dashed lines) if the pase is AF are obtained. The inset shows the behavior of 
the typical convergence time scale as a function of the distance to the critical 
point 5J±_ — J± — J1. The curves are power law fits (see text). 



\m{p + 1) — m{p)\ which is exponentially vanishing 

V{p) cx exp(-i) for t = ^ >> 1, (44) 

ST 

as we can see in Fig. |H1 It is very important to note that the convergence 
of the MF procedure is universal in the sense that the choice of the starting 
value to(0) is not crucial. 

The time scale is J^-dependant and diverges as a power law ^ l^Jl^'^ 
when approaching the critical line where \6J\ = | Jj_ — Jj^|. Our datas suggest 
/X ~ 1.06 if Jl > Jj_ and n ~ 0.95 if > for L = 12 (see inset of Fig.|HIl. 
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4 Lanczos algorithm for non uniform 
systems : Application to doped SP chains 

Doping a SP system with non-magnetic impurities leads to very surprising 
new features. For example in Cui_2;Ma;Ge03 (M==Zn or Mg), the discovery of 
coexistence between dimerization and AF long range order at small impurity 
concentration has motivated extented experimental j39j and theoretical 351 
I4UI I41L Pm 1431 144j investigations. In the following we report numerical studies 
of models for doped coupled spin chains. For sake of completness we also 
include a four-spin coupling which originates from cyclic exchange 0^1 • 

4.1 Doped coupled frustrated spin-^ chain with four-spin 
exchange 

As for the transverse coupling J± in Eas. H40l41() . we also apply the MF treat- 
ment to the added 4-spin coupling J4{Si^a ■ Si+i^a){Si+i,a+i ■ Si^a+i)- This 
leads to a self-consistent modulation of the NN couplings 

Hcs{a, J_l, J4) = X/t^-*- + SJi^a)Si,a ■ Si+l,a 
i,a 

+aSt^a ■ Si+2,a + hi^aSla\ + coustaut , (45) 
with hi^a given by Ea. (|41(l and 

5Ji,a — Ji{{Si^a+l ■ S'i+l.a+l) + ' Si+i^a-l)}- (46) 

Such a modulation produced by the J4 term stabilizes the SP phase and 
raises the transition line in Figd (for more details, see 02]). Another inter- 
esting feature of this model is its direct link with the magneto-elastic model 
considered in |43j where the elastic coupling K plays a role very similar to 
that of I/J4. In the following the parameters a, J± and J4 are set in order 
to constrain the system to be in a SP state in the absence of dopants. 

A dopant is described as an inert site i.e. all couplings to this site will 
be set to zero. Contrary to what we have seen previously, the use of the 
translation invariance is now forbidden by the presence of a single defects or 
by randomly located defects (see Figl^)). The maximal size accessible with 
a Lanczos procedure is then reduced because of this lack of symmetry and 
also because of the repeated iterative MF procedure. Indeed, the problem 
can not be reduced to a single chain model and hence the M non- equivalent 
chains have to be diagonalized independently. Following the method used 
in Ref . |44| . the MF equations are solved self-consistently on finite L x M 
clusters. Therefore, in the doped case, the time scale of the MF convergence 
^1- for the single chain problem is typically multiplied by M. 
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Fig. 9. Schematic picture of the coupled chains model with nearest neighbor, next- 
nearest neighbor, inter-chain and 4-spin couplings Ji, J2 = aJi, J±, and J4. Full 
(resp. open) circles stand for spin-i sites (resp. non-magnetic dopants). 

4.2 Confinement 

Replacing a single spin-i in a spontaneously dimerized (isolated) spin chain 
by a non magnetic dopant (described as an inert site) liberates a free spin 
i, named a soliton, which does not bind to the dopant (35|. The soliton can 
be depicted as a single unpaired spin (domain) separating two dimcr con- 
figurations The physical picture is completely different when a static 
bond dimerisation exists and produces an attractive potential between the 
soliton and the dopant [SEIEIII and consequently leads, under doping, to the 
formation of local magnetic moments [Hfil I42j as well as a rapid suppression 
of the spin gap 01]. However, a coupling to a purely one-dimensional (ID) 
adiabatic lattice @3I does not produce confinement in contrast to more re- 
alistic models including an elastic inter-chain coupling (to mimic 2D or 3D 
lattices) gSESj. 

Here, we re-examine the confinement problem in the context of the previous 
model including interchain magnetic coupling. 

a) Different kinds of dimer orders 

Let us return to model (|45|l . For J4 — 0, the MF treatment of the transverse 
magnetic coupling J_l does not break the degeneracy of the ground state : each 
chain displays a 2-fold degenerate ground state (dimers can stand either on 
even or odd bonds) independently from the other ones. The situation changes 
radically when J4 7^ because the degeneracy is reduced to 2. Indeed, each 
chain displays the same dimerized pattern if J4 < (columnar dimer order) 
whereas the dimer order is staggered in the transverse direction if J4 > 0, as 
we can see in FigllOl Consequently, the soliton remains deconfined 03 when 
J4 = as we can observe in Fiers lllll^ On the other hand, if J4 ^ the bulk 
dimerization constrains the soliton to lie in the vicinity of the impurity (see 
Figini)- 
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Fig. 10. Energy difference Zistg-cim between ground states witli staggered and 
collumnar dimer orders plotted versus J4 for model I45II (without dopant) at a = 
0.5, Ji = 0.1 and L = 12. 












(a) J4 = 



(b) J4 < 



(c) J4 > 



Fig. 11. Schematic picture of the soliton confinement mechanism induced by the 
coupling J4 of the model 1451 . The non magnetic dopant is represented by an open 
circle and the large black bonds stand for stronger dimer bonds. The black arrow 
represents the soliton, released by the impurity, which is deconfined if J4 = (a) 
whereas it is linked to the dopant if J4 7^ 0. We can see that this binding is imposed 
by the bulk dimerization which is columnar if J4 < (b) or staggered if J4 > (c) 



h ) Enhancement of the magnetization near a dopant 

Under doping, the system becoming inhomogeneous, we define a local mean 
staggered magnetization 



(47) 



which has been calculated for a single dopant in a system of size L x M = 
16 X 8. It is plotted for different values of the four-spin coupling in Fig. [T^ 
where the confinement mechanism can clearly be observed. Note that the 
inter-chain coupling induces a "polarization cloud" with strong AF correla- 
tions in the neighbor chains of the doped one. 



Simulatic 
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• J4 = 




Fig. 12. Local magnetization Mf^^ for L x M=16 x 8 coupled chains with one 
dopant D (shown by arrow) located at a = 1, i = 16 in the dimerised phase {a = 0.5, 
Jj_ ~ 0.1). Circles correspond to J4 = (shown up to the third neighbor chain of 
the doped one) and squares (crosses) to J4 = 0.01 (J4 = 0.08). The coupling J2 
across the dopant has been set to for convenience (Reprinted from Ref.|47p. 

c) Confinement length 

In order to measure the strength of confinement, a confinement length can 
be defined as 

^ Ml' 

In the absence of confinement, the solitonic cloud is located at the center 
of the doped chain : ^|| = L/2. Otherwise, ^|| converges to a finite value 
when L 00. On Fig ll3l the confinement lenght is plotted versus J4 for 
2 different system sizes at a = 0.5 and J± ~ 0.1. The finite size effects 
decrease for increasing J4. Note that ^||(J4) 7^ C||(^>^4) and a power law ^U] 
with different exponents rj is expected when J4 ^ 0. A fit gives 77 ~ 0.33 if 
J4 < and 77 ~ 0.50 for J4 > fFig ll3|l . This asymmetry can be understood 
from opposite renormalisations of Ji for different signs of J4. Indeed, if J4 < 
then SJi^a > and the nearest neighbor MF exchange becomes larger than 
the bare one. Opposite effects are induced by J4 > 0. 

4.3 Effective interaction 

We now turn to the investigation of the effective interaction between dopants. 
A system of coupled chains with two dopants is considered here (see FiglHll. 
Each impurity releases an effective spin i, localized at a distance ~ ^|| from 
it due to the confining potentiel set by J4. We define an effective pairwise 
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Fig. 13. ED data of the soliton average position vs J4 calculated for a = 0.5 and 
= 0.1. Different symbols are used for L x M = 12 x 6 and 16 x 8 clusters. The 
long-dashed line is a power-law fit (see text). Inset shows the magnetization profile 
in the doped (a = 1) chain at J4 = 0.08, ie ~ 2.5 (Reprinted from Ref.|47p. 

interaction J'^^ as the energy difference of the S' = 1 and the 5* = ground 
states. When J'^^ = E{S = 1) — E{S = 0) is positive (negative) the spin 
interaction is AF (ferromagnetic). Let us first consider the case of two dopants 
in the same chain, (i) When the two vacancies are on the same sub-lattice the 
moments experience a very small ferromagnetic < as seen in Fig. 1141 
with Aa = so that the two effective spins ^ are almost free, (ii) When the 
two vacancies sit on different sub-lattices, Ai is odd and the effective coupling 
is AF with a magnitude close to the singlet-triplet gap. Fig. E| with Aa — 
shows that the decay of J'^^ with distance is in fact very slow for such a 
configuration. Physically, this result shows that a soliton and an anti-soliton 
on the same chain and different sublattices tend to recombine. 

The behavior of the pairwise interaction of two dopants located on differ- 
ent chains {Aa = 1, 2, 3, 4) is shown on Fig. E|for Aa = 1, 2, 3, 4 for J4 > 0. 
When dopants are on opposite sub-lattices the effective interaction is antifer- 
romagnetic. At small dopant separation J'^^{Ai) increases with the dopant 
separation as the overlap between the two AF clouds increases until Ai ~ 2^. 
For larger separation, J°^{Ai) decays rapidly. Note that the released spin-i 
solitons bind on the opposite right and left sides of the dopants as imposed 
by the the bulk dimerisation flS]. If dopants are on the same sub- lattice, 
solitons are located on the same side of the dopants (SHI and the effective 
exchange J'^^{Ai) is ferromagnetic and decays rapidly to become negligible 
when Ai > 2^. The key feature here is the fact that the effective pairwise 
interaction is not frustrating (because of its sign alternation with distance) 
although frustration is present in the microscopic underlying model. AF or- 
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0.01 



Fig. 14. Magnitude of the effective magnetic coupling between two impurities 
located either on the same chain (a) or on different ones (b-c-d) vs the dopant 
separation Ai in a system of size L x M — 16 x 8 with a = 0.5, J± — 0.1, and 
J4 = 0.08. Closed (resp. open) symbols correspond to AF (F) interactions. 



dering is then expected (at T = 0) as seen for a related system of coupled 
Spin-Peierls chains |44| . 



5 Conclusion 

The coexistence between AF order and SP dimer order under doping SP ma- 
terials with non magnetic impurities |39j is one of the most surprising phe- 
nomenon in the field of quantum magnetism. Starting with the non frustating 
interaction between two solitonic clouds calculated above, we can construct 
an effective model of long range interacting spins ^ , randomly diluted on a 
square lattice. We have implemented a Quantum Monte Carlo (QMC) algo- 
rithm using the Stochastic Series Expansion (SSE) method |H3 in order to 
study long distance interacting spin-i models. The mechanism of AF order- 
ing has been studied at very low temperature and dopant concentration with 
very large system sizes, up to 96 x 96 |52| . 
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